#### Replication code for "The Black Market Blues: The Political Costs of Illicit Currency Markets"

rm(list = ls())
setwd("~/Dropbox/Argentina Black Market/JOP Final Submission")

# Packages
library(list)
library(xtable)
library(ggplot2)
library(gridExtra)


load('Black_Market.RData')


# Table E1 Appendix

desc<-rbind(c(sum(!is.na(df$pj)),mean(df$pj,na.rm=T),sd(df$pj,na.rm=T),min(df$pj,na.rm=T),max(df$pj,na.rm=T)),
            c(sum(!is.na(df$cambiemos)),mean(df$cambiemos,na.rm=T),sd(df$cambiemos,na.rm=T),
              min(df$cambiemos,na.rm=T),max(df$cambiemos,na.rm=T)),
            c(sum(!is.na(df$ideology)),mean(df$ideology,na.rm=T),sd(df$ideology,na.rm=T),
              min(df$ideology,na.rm=T),max(df$ideology,na.rm=T)),
            c(sum(!is.na(df$middleclass)),mean(df$middleclass,na.rm=T),sd(df$middleclass,na.rm=T),
              min(df$middleclass,na.rm=T),max(df$middleclass,na.rm=T)),
            c(sum(!is.na(df$upperclass)),mean(df$upperclass,na.rm=T),sd(df$upperclass,na.rm=T),
              min(df$upperclass,na.rm=T),max(df$upperclass,na.rm=T)),
            c(sum(!is.na(df$college)),mean(df$college,na.rm=T),sd(df$college,na.rm=T),
              min(df$college,na.rm=T),max(df$college,na.rm=T)),
            c(sum(!is.na(df$female)),mean(df$female,na.rm=T),sd(df$female,na.rm=T),
              min(df$female,na.rm=T),max(df$female,na.rm=T)),
            c(sum(!is.na(df$age)),mean(df$age,na.rm=T),sd(df$age,na.rm=T),
              min(df$age,na.rm=T),max(df$age,na.rm=T)),
            c(sum(!is.na(df$amba)),mean(df$amba,na.rm=T),sd(df$amba,na.rm=T),
              min(df$amba,na.rm=T),max(df$amba,na.rm=T)),
            c(sum(!is.na(df$sciolivote)),mean(df$sciolivote,na.rm=T),sd(df$sciolivote,na.rm=T),
              min(df$sciolivote,na.rm=T),max(df$sciolivote,na.rm=T)),
            c(sum(!is.na(df$macri)),mean(df$macri,na.rm=T),sd(df$macri,na.rm=T),
              min(df$macri,na.rm=T),max(df$macri,na.rm=T)),
            c(sum(!is.na(df$cepo)),mean(df$cepo,na.rm=T),sd(df$cepo,na.rm=T),
              min(df$cepo,na.rm=T),max(df$cepo,na.rm=T)),
            c(sum(!is.na(df$know_cepo)),mean(df$know_cepo,na.rm=T),sd(df$know_cepo,na.rm=T),
              min(df$know_cepo,na.rm=T),max(df$know_cepo,na.rm=T)),
            c(sum(!is.na(df$know_cepo2)),mean(df$know_cepo2,na.rm=T),sd(df$know_cepo2,na.rm=T),
              min(df$know_cepo2,na.rm=T),max(df$know_cepo2,na.rm=T)),
            c(sum(!is.na(df$know_ypf)),mean(df$know_ypf,na.rm=T),sd(df$know_ypf,na.rm=T),
              min(df$know_ypf,na.rm=T),max(df$know_ypf,na.rm=T)),
            c(sum(!is.na(df$know_auh)),mean(df$know_auh,na.rm=T),sd(df$know_auh,na.rm=T),
              min(df$know_auh,na.rm=T),max(df$know_auh,na.rm=T)),
            c(sum(!is.na(df$know_gay)),mean(df$know_gay,na.rm=T),sd(df$know_gay,na.rm=T),
              min(df$know_gay,na.rm=T),max(df$know_gay,na.rm=T)))



rows<-c('Peronist Party ID','Cambiemos Party ID','Ideology','Middle class',
        'Upper class','College','Female','Age','Amba','Pr(Scioli Vote)',
        'Approval: Macri','Approval: Cepo', 'Knowledge: cepo adoption',
        'Knowledge: cepo removal','Knowledge: YPF nationalization',
        'Knowledge: CCTs', 'Knowledge: Gay marriage law' )

desc<-data.frame(cbind(rows,desc))


colnames(desc)<-c('Variable','Observations','Mean','Std. dev.', 'Min','Max')
desc<-xtable(desc,digits=2,type='text')
print(desc, include.rownames = F,sanitize.text.function=function(x){x})


########## Appendix Table E2

# Upper panel: 
desc2<-data.frame(0:5,cbind(c(with(df,table(list4)),''),
                           as.numeric(c(100*with(df,prop.table(table(list4))),'')),
                           with(df,table(list5)),100*with(df,prop.table(table(list5)))))
colnames(desc2)<-c('Number of activities','Control group','Treatment group')

print(desc2)
# Lower panel: Means and N for treatment and control groups were added afterwards

mean(df$list_outcome[df$listgroup==0],na.rm=T)  # Mean control
mean(df$list_outcome[df$listgroup==1],na.rm=T)  # Mean treatment
table(df$listgroup) # N for each group



# Testing for Design effects (Reported in Appendix)
test <- ict.test(y = df$list_outcome[!is.na(df$list_outcome)], treat = 
                   df$listgroup[!is.na(df$list_outcome)], J = 4, alpha = 0.05, gms = FALSE)
test

########## Core Full MLE Model: Scioli vote as outcome


scioli<- ictreg.joint(formula = list_outcome ~ pj+cambiemos+ideology+middleclass+upperclass+college+female+age+
                        amba, data = df,
                      treat = "listgroup", outcome = "sciolivote", J = 4, constrained = TRUE,
                      outcome.reg = "logistic", maxIter = 1000)



################################# Figure 2: Who buys dollars in the black market? ###########################################



# Predicted probabilities of buying dollars in black market by college degree, class, and residence



edul.no.ba<-predict.ictreg.joint(scioli,newdata=subset(df,college==0 & upperclass==0 & middleclass ==0 & amba==1),
                                 se.fit = TRUE, interval = "confidence", level = 0.9,
                                 avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)

edul.yes.ba<-predict.ictreg.joint(scioli,newdata=subset(df,college==1 & upperclass==0 & middleclass ==0& amba==1),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = TRUE)


edum.no.ba<-predict.ictreg.joint(scioli,newdata=subset(df,college==0 &  middleclass ==1& amba==1),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)


edum.yes.ba<-predict.ictreg.joint(scioli,newdata=subset(df,college==1 &  middleclass ==1& amba==1),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = TRUE)

eduu.no.ba<-predict.ictreg.joint(scioli,newdata=subset(df,college==0 &  upperclass ==1& amba==1),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)


eduu.yes.ba<-predict.ictreg.joint(scioli,newdata=subset(df,college==1 &  upperclass ==1& amba==1),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = TRUE)






edul.no.in<-predict.ictreg.joint(scioli,newdata=subset(df,college==0 & upperclass==0 & middleclass ==0 & amba==0),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)

edul.yes.in<-predict.ictreg.joint(scioli,newdata=subset(df,college==1 & upperclass==0 & middleclass ==0& amba==0),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = TRUE)


edum.no.in<-predict.ictreg.joint(scioli,newdata=subset(df,college==0 &  middleclass ==1& amba==0),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)


edum.yes.in<-predict.ictreg.joint(scioli,newdata=subset(df,college==1 &  middleclass ==1& amba==0),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = TRUE)

eduu.no.in<-predict.ictreg.joint(scioli,newdata=subset(df,college==0 &  upperclass ==1& amba==0),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)


eduu.yes.in<-predict.ictreg.joint(scioli,newdata=subset(df,college==1 &  upperclass ==1& amba==0),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = TRUE)

plot.who<-rbind(edul.no.ba$fitsens,edul.yes.ba$fitsens,edum.no.ba$fitsens,edum.yes.ba$fitsens,eduu.no.ba$fitsens,eduu.yes.ba$fitsens,
                edul.no.in$fitsens,edul.yes.in$fitsens,edum.no.in$fitsens,edum.yes.in$fitsens,eduu.no.in$fitsens,eduu.yes.in$fitsens)



plot.who$edu<-rep(c(c('No','Yes')),6)
plot.who$Class<-c(rep('Lower',2),rep('Middle',2),rep('Upper',2),rep('Lower',2),rep('Middle',2),rep('Upper',2))

plot.who$region<-c(rep('AMBA',6),rep('Non-AMBA',6))


pd <- position_dodge(.5)


ggplot(plot.who,aes(x=Class,y=fit,colour=edu,shape=edu))+facet_wrap(~region)+
  geom_hline(aes(yintercept=0), colour="#999999", linetype="dashed") +
  geom_pointrange( data =plot.who, stat = "identity", size =.7,position=pd,
                   mapping = aes(ymin= lwr,ymax=upr))+
  ylab("Proportion buying dollars in black market") +xlab("Class self-perception") +
  theme_bw() +
  #theme(aspect.ratio=.8)+
  theme(legend.position=c(.1,.9))+
  theme(axis.text.x  = element_text( size=12),axis.title.x = element_text(size=14)) +
  theme(axis.text.y  = element_text( size=12),axis.title.y = element_text(size=14))+
  theme(legend.background = element_rect( fill="white", size=.5, linetype="dotted"))+
  theme(legend.title = element_text(colour = 'black', size = 10))+
  theme(legend.text = element_text(colour = 'black', size = 8,vjust=0, hjust = 0))+
  theme(legend.key.width = unit(.5, "cm"))+
  theme(legend.key.size = unit(.5, "cm"))+
  theme(plot.margin = unit(c(.5, .5,.5, .5), "cm"))+
  
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank())+theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+
  scale_y_continuous(breaks=seq(0,.75,.25),labels=seq(0,.75,.25),lim=c(0,.75))+
  guides(shape=guide_legend(title="College Degree"),colour=guide_legend(title="College Degree"))



# Appendix Figure E1: Comparing predictions of who buy dollars from models with and without
# correction for floor and ceiling effects



# MLE model of who buys dollars with correction for floor/ceiling effects


m2 <- ictreg(list_outcome ~ pj+cambiemos+ideology+middleclass+upperclass+college+female+age+amba+
               amba, treat = "listgroup",
             J = 4, data = df, method = "ml", fit.start = "glm", 
             floor = TRUE, ceiling = TRUE, 
             floor.fit = "bayesglm", ceiling.fit = "bayesglm",
             floor.formula = ~pj+cambiemos+ideology+middleclass+upperclass+college+female+age+amba,
             ceiling.formula = ~ pj+cambiemos+ideology+middleclass+upperclass+college+female+age+amba)


edul.no.ba<-predict.ictreg(m2,newdata=subset(df,college==0 & upperclass==0 & middleclass ==0 & amba==1),
                           se.fit = TRUE, interval = "confidence", level = 0.95,
                           avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                           predict.sensitive = TRUE)

edul.yes.ba<-predict.ictreg(m2,newdata=subset(df,college==1 & upperclass==0 & middleclass ==0& amba==1),
                            se.fit = TRUE, interval = "confidence", level = 0.95,
                            avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                            predict.sensitive = TRUE)


edum.no.ba<-predict.ictreg(m2,newdata=subset(df,college==0 &  middleclass ==1& amba==1),
                           se.fit = TRUE, interval = "confidence", level = 0.95,
                           avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                           predict.sensitive = TRUE)


edum.yes.ba<-predict.ictreg(m2,newdata=subset(df,college==1 &  middleclass ==1& amba==1),
                            se.fit = TRUE, interval = "confidence", level = 0.95,
                            avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                            predict.sensitive = TRUE)

eduu.no.ba<-predict.ictreg(m2,newdata=subset(df,college==0 &  upperclass ==1& amba==1),
                           se.fit = TRUE, interval = "confidence", level = 0.95,
                           avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                           predict.sensitive = TRUE)


eduu.yes.ba<-predict.ictreg(m2,newdata=subset(df,college==1 &  upperclass ==1& amba==1),
                            se.fit = TRUE, interval = "confidence", level = 0.95,
                            avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                            predict.sensitive = TRUE)






edul.no.in<-predict.ictreg(m2,newdata=subset(df,college==0 & upperclass==0 & middleclass ==0 & amba==0),
                           se.fit = TRUE, interval = "confidence", level = 0.95,
                           avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                           predict.sensitive = TRUE)

edul.yes.in<-predict.ictreg(m2,newdata=subset(df,college==1 & upperclass==0 & middleclass ==0& amba==0),
                            se.fit = TRUE, interval = "confidence", level = 0.95,
                            avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                            predict.sensitive = TRUE)


edum.no.in<-predict.ictreg(m2,newdata=subset(df,college==0 &  middleclass ==1& amba==0),
                           se.fit = TRUE, interval = "confidence", level = 0.95,
                           avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                           predict.sensitive = TRUE)


edum.yes.in<-predict.ictreg(m2,newdata=subset(df,college==1 &  middleclass ==1& amba==0),
                            se.fit = TRUE, interval = "confidence", level = 0.95,
                            avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                            predict.sensitive = TRUE)

eduu.no.in<-predict.ictreg(m2,newdata=subset(df,college==0 &  upperclass ==1& amba==0),
                           se.fit = TRUE, interval = "confidence", level = 0.95,
                           avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                           predict.sensitive = TRUE)


eduu.yes.in<-predict.ictreg(m2,newdata=subset(df,college==1 &  upperclass ==1& amba==0),
                            se.fit = TRUE, interval = "confidence", level = 0.95,
                            avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                            predict.sensitive = TRUE)

data.ict2<-rbind(edul.no.ba$fit,edul.yes.ba$fit,edum.no.ba$fit,edum.yes.ba$fit,eduu.no.ba$fit,eduu.yes.ba$fit,
                edul.no.in$fit,edul.yes.in$fit,edum.no.in$fit,edum.yes.in$fit,eduu.no.in$fit,eduu.yes.in$fit)



data.ict2$edu<-rep(c(c('No','Yes')),6)
data.ict2$Class<-c(rep('Lower',2),rep('Middle',2),rep('Upper',2),rep('Lower',2),rep('Middle',2),rep('Upper',2))

data.ict2$region<-c(rep('AMBA',6),rep('Non-AMBA',6))



edul.no.ba<-predict.ictreg.joint(scioli,newdata=subset(df,college==0 & upperclass==0 & middleclass ==0 & amba==1),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)

edul.yes.ba<-predict.ictreg.joint(scioli,newdata=subset(df,college==1 & upperclass==0 & middleclass ==0& amba==1),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = TRUE)


edum.no.ba<-predict.ictreg.joint(scioli,newdata=subset(df,college==0 &  middleclass ==1& amba==1),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)


edum.yes.ba<-predict.ictreg.joint(scioli,newdata=subset(df,college==1 &  middleclass ==1& amba==1),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = TRUE)

eduu.no.ba<-predict.ictreg.joint(scioli,newdata=subset(df,college==0 &  upperclass ==1& amba==1),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)


eduu.yes.ba<-predict.ictreg.joint(scioli,newdata=subset(df,college==1 &  upperclass ==1& amba==1),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = TRUE)






edul.no.in<-predict.ictreg.joint(scioli,newdata=subset(df,college==0 & upperclass==0 & middleclass ==0 & amba==0),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)

edul.yes.in<-predict.ictreg.joint(scioli,newdata=subset(df,college==1 & upperclass==0 & middleclass ==0& amba==0),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = TRUE)


edum.no.in<-predict.ictreg.joint(scioli,newdata=subset(df,college==0 &  middleclass ==1& amba==0),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)


edum.yes.in<-predict.ictreg.joint(scioli,newdata=subset(df,college==1 &  middleclass ==1& amba==0),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = TRUE)

eduu.no.in<-predict.ictreg.joint(scioli,newdata=subset(df,college==0 &  upperclass ==1& amba==0),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)




eduu.yes.in<-predict.ictreg.joint(scioli,newdata=subset(df,college==1 &  upperclass ==1& amba==0),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = TRUE)

data.ict<-rbind(edul.no.ba$fitsens,edul.yes.ba$fitsens,edum.no.ba$fitsens,edum.yes.ba$fitsens,eduu.no.ba$fitsens,eduu.yes.ba$fitsens,
                edul.no.in$fitsens,edul.yes.in$fitsens,edum.no.in$fitsens,edum.yes.in$fitsens,eduu.no.in$fitsens,eduu.yes.in$fitsens)



data.ict$edu<-rep(c(c('No','Yes')),6)
data.ict$Class<-c(rep('Lower',2),rep('Middle',2),rep('Upper',2),rep('Lower',2),rep('Middle',2),rep('Upper',2))

data.ict$region<-c(rep('AMBA',6),rep('Non-AMBA',6))

data.plot<-rbind(data.ict2,data.ict)
data.plot$ict<-c(rep('Floor and Ceiling Effects',12),rep('No correction',12))


# Appendix Figure E1
ggplot(data=data.plot,aes(x=Class,y=fit,colour=edu,shape=ict,linetype=edu))+facet_wrap(~region)+
  geom_hline(aes(yintercept=0), colour="#999999", linetype="dashed") +
  geom_pointrange( data =data.plot, stat = "identity", size =.5,position=pd,
                   mapping = aes(ymin= lwr,ymax=upr))+
  ylab("Proportion buying dollars in black market") +xlab("Class self-perception") +
  theme_bw() +
  theme(axis.text.x  = element_text( size=12),axis.title.x = element_text(size=14)) +
  theme(axis.text.y  = element_text( size=12),axis.title.y = element_text(size=14))+
  theme(legend.background = element_rect( fill="white", size=.5, linetype="dotted"))+
  theme(legend.title = element_text(colour = 'black', size = 6))+
  theme(legend.text = element_text(colour = 'black', size = 6,vjust=0, hjust = 0))+
  theme(legend.key.width = unit(.5, "cm"))+
  theme(legend.key.size = unit(.5, "cm"))+
  theme(plot.margin = unit(c(.5, .5,.5, .5), "cm"))+
  
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank())+theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+
  guides(shape=guide_legend(title="Correction for lyers?",ncol=2),linetype=guide_legend(title="College Degree",ncol=2),
         colour=guide_legend(title="College Degree",ncol=2))+
  theme(legend.position = c(.65,.85))




### Appendix Figure E2. List experiment. Balance tests.


edu<-with(df,t.test(college~listgroup))
age<-with(df,t.test(age~listgroup))
female<-with(df,t.test(female~listgroup))
middle<-with(df,t.test(middleclass~listgroup))
upper<-with(df,t.test(upperclass~listgroup))
pj<-with(df,t.test(pj~listgroup))
cambiemos<-with(df,t.test(cambiemos~listgroup))
ideology<-with(df,t.test(ideology~listgroup))
sciolivote<-with(df,t.test(sciolivote~listgroup))
macri<-with(df,t.test(macri~listgroup))
cepo<-with(df,t.test(cepo~listgroup))
cepo2<-with(df,t.test(know_cepo~listgroup))
cepo3<-with(df,t.test(know_cepo2~listgroup))


balance <- rbind(cbind('College degree',sum(!is.na(df$college)),edu$estimate[1],edu$estimate[2],
                       edu$p.value,with(df,wilcox.test(college~listgroup))[3]),
                 cbind('Age',sum(!is.na(df$age)),age$estimate[1],age$estimate[2],age$p.value,with(df,wilcox.test(age~listgroup))[3]),
                 cbind('Female',sum(!is.na(df$female)),female$estimate[1],edu$estimate[2],female$p.value,with(df,wilcox.test(female~listgroup))[3]),
                 cbind('Middle class',sum(!is.na(df$middleclass)),middle$estimate[1],middle$estimate[2],middle$p.value,with(df,wilcox.test(middleclass~listgroup))[3]),
                 cbind('Upper class',sum(!is.na(df$upperclass)),upper$estimate[1],upper$estimate[2],upper$p.value,with(df,wilcox.test(upperclass~listgroup))[3]),
                 cbind('PJ identifier',sum(!is.na(df$pj)),pj$estimate[1],pj$estimate[2],pj$p.value,with(df,wilcox.test(pj~listgroup))[3]),
                 cbind('Cambiemos identifier',sum(!is.na(df$cambiemos)),cambiemos$estimate[1],cambiemos$estimate[2],cambiemos$p.value,
                       with(df,wilcox.test(cambiemos~listgroup))[3]),
                 cbind('Ideology',sum(!is.na(df$ideology)),ideology$estimate[1],ideology$estimate[2],ideology$p.value,with(df,wilcox.test(ideology~listgroup))[3]),
                 cbind('Scioli Vote',sum(!is.na(df$sciolivote)),sciolivote$estimate[1],sciolivote$estimate[2],sciolivote$p.value,with(df,wilcox.test(sciolivote~listgroup))[3]),
                 cbind('Approval: President Macri',sum(!is.na(df$macri)),macri$estimate[1],macri$estimate[2],macri$p.value,
                       with(df,wilcox.test(macri~listgroup))[3]),
                 cbind('Approval: Cepo removal',sum(!is.na(df$cepo)),cepo$estimate[1],cepo$estimate[2],cepo$p.value,with(df,wilcox.test(cepo~listgroup))[3]),
                 cbind('Knowledge: Cepo adoption',sum(!is.na(df$know_cepo)),cepo2$estimate[1],cepo2$estimate[2],cepo2$p.value,
                       with(df,wilcox.test(know_cepo~listgroup))[3]),
                 cbind('Knowledge: Cepo removal',sum(!is.na(df$know_cepo2)),cepo3$estimate[1],cepo3$estimate[2],cepo3$p.value,with(df,wilcox.test(know_cepo2~listgroup))[3]))


balance <- matrix(data = balance, nrow = 13, ncol = 6, byrow = FALSE, dimnames = NULL)
balance[,3:6]<-round(as.numeric(balance[,3:6]),2)
results <- balance


pdf(file="Balance.pdf",family="Palatino", width=7,height=6.5)

textsize=0.8
parcex=0.9
at1=-0.35
at2=-0.15
at3=-0.9
xlim1=-0.85

xlim = c(xlim1,1); pchset = c(21,24,22,23); pchcolset = c("blue","yellow","red","darkgreen")

par(cex=parcex, mai = c(0.6, 0.35, 1.1, 0.35))

ny = nrow(results)

plot(x=NULL,axes=F, xlim=xlim, ylim=c(1,ny),xlab="",ylab="")

abline(v=c(0,0.05,0.1),lty=c(1,2,2), lwd=c(1,1,1))
axis(side=1,at=c(0,0.05,0.1,1),tick=TRUE, las=2, cex.axis=0.7)

axis(side=3,at=at1,labels="Mean\nTreated",tick=FALSE, padj=0.5,cex.axis=textsize)
axis(side=3,at=at2,labels="Mean\nUntreated",tick=FALSE, padj=0.5,cex.axis=textsize)
axis(side=3,at=0.5,labels="P-values",tick=FALSE, padj=0.5,cex.axis=textsize)

points(results[,5],ny:1, pch = 21, col = "blue", bg = "blue", cex=.7)
points(results[,6],ny:1, pch = 24, col = "red", bg = "red", cex=.7)

for(i in 1:ny) {
  text(at3,ny-i+1,results[i,1],adj = 0,cex=textsize) # variable name
  text(at1,ny-i+1,results[i,3], cex=textsize)        # treatment mean
  text(at2,ny-i+1,results[i,4], cex=textsize)        # control mean
}



legend(x=.75, y=2.2, c("t-test",'Rank-sum test'), pch=c(21,24), pt.bg = c("blue",'red'), cex=0.8, bty="n")

dev.off()




#### Appendix Figure E3. Marginal effect of covariates on probability of buying black market dollars

ci<-.9
low<-predict.ictreg.joint(scioli,newdata=subset(df,middleclass==0 & upperclass==0  ),
                          se.fit = TRUE, interval = "confidence", level = ci,
                          avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                          predict.sensitive = TRUE)



mid<-predict.ictreg.joint(scioli,newdata=subset(df,middleclass==1),
                          se.fit = TRUE, interval = "confidence", level = ci,
                          avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                          predict.sensitive = TRUE)

up<-predict.ictreg.joint(scioli,newdata=subset(df,upperclass==1),
                         se.fit = TRUE, interval = "confidence", level = ci,
                         avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                         predict.sensitive = TRUE)


predict.mid <- mid$draws.mean.sens - low$draws.mean.sens 
mean <- mean(predict.mid, na.rm = TRUE)
ci.low.95 <- quantile(predict.mid, probs = .025, na.rm = TRUE)
ci.upper.95 <- quantile(predict.mid, probs = .975, na.rm = TRUE)

mid.95 <- cbind(ci.low.95, mean, ci.upper.95)

ci.low.90 <- quantile(predict.mid, probs = .05, na.rm = TRUE)
ci.upper.90 <- quantile(predict.mid, probs = .95, na.rm = TRUE)
mid.90 <- cbind(ci.low.90, mean, ci.upper.90)




predict.up <- up$draws.mean.sens - low$draws.mean.sens 
mean <- mean(predict.up, na.rm = TRUE)
ci.low.95 <- quantile(predict.up, probs = .025, na.rm = TRUE)
ci.upper.95 <- quantile(predict.up, probs = .975, na.rm = TRUE)

up.95 <- cbind(ci.low.95, mean, ci.upper.95)
up

ci.low.90 <- quantile(predict.up, probs = .05, na.rm = TRUE)
ci.upper.90 <- quantile(predict.up, probs = .95, na.rm = TRUE)
up.90 <- cbind(ci.low.90, mean, ci.upper.90)





int<-predict.ictreg.joint(scioli,newdata=subset(df,amba==0  ),
                          se.fit = TRUE, interval = "confidence", level = ci,
                          avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                          predict.sensitive = TRUE)

ba<-predict.ictreg.joint(scioli,newdata=subset(df,amba==1),
                         se.fit = TRUE, interval = "confidence", level = ci,
                         avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                         predict.sensitive = TRUE)


predict.ba <- ba$draws.mean.sens - int$draws.mean.sens 
mean <- mean(predict.ba, na.rm = TRUE)
ci.low.95 <- quantile(predict.ba, probs = .025, na.rm = TRUE)
ci.upper.95 <- quantile(predict.ba, probs = .975, na.rm = TRUE)

ba.95 <- cbind(ci.low.95, mean, ci.upper.95)


ci.low.90 <- quantile(predict.ba, probs = .05, na.rm = TRUE)
ci.upper.90 <- quantile(predict.ba, probs = .95, na.rm = TRUE)

ba.90 <- cbind(ci.low.90, mean, ci.upper.90)






col<-predict.ictreg.joint(scioli,newdata=subset(df,college==1  ),
                          se.fit = TRUE, interval = "confidence", level = ci,
                          avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                          predict.sensitive = TRUE)

no<-predict.ictreg.joint(scioli,newdata=subset(df,college==0),
                         se.fit = TRUE, interval = "confidence", level = ci,
                         avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                         predict.sensitive = TRUE)


predict.col <- col$draws.mean.sens - no$draws.mean.sens 
mean <- mean(predict.col, na.rm = TRUE)
ci.low.95 <- quantile(predict.col, probs = .025, na.rm = TRUE)
ci.upper.95 <- quantile(predict.col, probs = .975, na.rm = TRUE)

col.95<- cbind(ci.low.95, mean, ci.upper.95)


ci.low.90 <- quantile(predict.col, probs = .05, na.rm = TRUE)
ci.upper.90 <- quantile(predict.col, probs = .95, na.rm = TRUE)

col.90<- cbind(ci.low.90, mean, ci.upper.90)
col




pj<-predict.ictreg.joint(scioli,newdata=subset(df,pj==1  ),
                         se.fit = TRUE, interval = "confidence", level = ci,
                         avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                         predict.sensitive = TRUE)

npj<-predict.ictreg.joint(scioli,newdata=subset(df,pj==0),
                          se.fit = TRUE, interval = "confidence", level = ci,
                          avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                          predict.sensitive = TRUE)


predict.pj <- pj$draws.mean.sens - npj$draws.mean.sens 
mean <- mean(predict.pj, na.rm = TRUE)
ci.low.95 <- quantile(predict.pj, probs = .025, na.rm = TRUE)
ci.upper.95 <- quantile(predict.pj, probs = .975, na.rm = TRUE)

pj.95 <- cbind(ci.low.95, mean, ci.upper.95)


ci.low.90 <- quantile(predict.pj, probs = .05, na.rm = TRUE)
ci.upper.90 <- quantile(predict.pj, probs = .95, na.rm = TRUE)

pj.90 <- cbind(ci.low.90, mean, ci.upper.90)




cambiemos<-predict.ictreg.joint(scioli,newdata=subset(df,cambiemos==1  ),
                                se.fit = TRUE, interval = "confidence", level = ci,
                                avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                predict.sensitive = TRUE)

ncambiemos<-predict.ictreg.joint(scioli,newdata=subset(df,cambiemos==0),
                                 se.fit = TRUE, interval = "confidence", level = ci,
                                 avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)


predict.cambiemos <- cambiemos$draws.mean.sens - ncambiemos$draws.mean.sens 
mean <- mean(predict.cambiemos, na.rm = TRUE)
ci.low.95 <- quantile(predict.cambiemos, probs = .025, na.rm = TRUE)
ci.upper.95 <- quantile(predict.cambiemos, probs = .975, na.rm = TRUE)

cambiemos.95 <- cbind(ci.low.95, mean, ci.upper.95)
cambiemos

ci.low.90 <- quantile(predict.cambiemos, probs = .05, na.rm = TRUE)
ci.upper.90 <- quantile(predict.cambiemos, probs = .95, na.rm = TRUE)

cambiemos.90 <- cbind(ci.low.90, mean, ci.upper.90)



female<-predict.ictreg.joint(scioli,newdata=subset(df,female==1  ),
                             se.fit = TRUE, interval = "confidence", level = ci,
                             avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                             predict.sensitive = TRUE)

nfemale<-predict.ictreg.joint(scioli,newdata=subset(df,female==0),
                              se.fit = TRUE, interval = "confidence", level = ci,
                              avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                              predict.sensitive = TRUE)


predict.female <- female$draws.mean.sens - nfemale$draws.mean.sens 
mean <- mean(predict.female, na.rm = TRUE)
ci.low.95 <- quantile(predict.female, probs = .025, na.rm = TRUE)
ci.upper.95 <- quantile(predict.female, probs = .975, na.rm = TRUE)

female.95 <- cbind(ci.low.95, mean, ci.upper.95)
female


ci.low.90 <- quantile(predict.female, probs = .05, na.rm = TRUE)
ci.upper.90 <- quantile(predict.female, probs = .95, na.rm = TRUE)

female.90 <- cbind(ci.low.90, mean, ci.upper.90)


age<-predict.ictreg.joint(scioli,newdata=subset(df,age==3  ),
                           se.fit = TRUE, interval = "confidence", level = ci,
                           avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                           predict.sensitive = TRUE)

age1<-predict.ictreg.joint(scioli,newdata=subset(df,age==1),
                           se.fit = TRUE, interval = "confidence", level = ci,
                           avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                           predict.sensitive = TRUE)


predict.age <- age$draws.mean.sens - age1$draws.mean.sens 
mean <- mean(predict.age, na.rm = TRUE)
ci.low.95 <- quantile(predict.age, probs = .025, na.rm = TRUE)
ci.upper.95 <- quantile(predict.age, probs = .975, na.rm = TRUE)

age.95 <- cbind(ci.low.95, mean, ci.upper.95)

ci.low.90 <- quantile(predict.age, probs = .05, na.rm = TRUE)
ci.upper.90 <- quantile(predict.age, probs = .95, na.rm = TRUE)

age.90 <- cbind(ci.low.90, mean, ci.upper.90)

ideology<-predict.ictreg.joint(scioli,newdata=subset(df,ideology==10  ),
                               se.fit = TRUE, interval = "confidence", level = ci,
                               avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                               predict.sensitive = TRUE)

nid<-predict.ictreg.joint(scioli,newdata=subset(df,ideology==0),
                          se.fit = TRUE, interval = "confidence", level = ci,
                          avg = TRUE, sensitive.value = "both", sensitive.diff = TRUE, return.draws = TRUE,
                          predict.sensitive = TRUE)


predict.ideology <- ideology$draws.mean.sens - nid$draws.mean.sens 
mean <- mean(predict.ideology, na.rm = TRUE)
ci.low.95 <- quantile(predict.ideology, probs = .025, na.rm = TRUE)
ci.upper.95 <- quantile(predict.ideology, probs = .975, na.rm = TRUE)

ideology.95 <- cbind(ci.low.95, mean, ci.upper.95)
ideology

ci.low.90 <- quantile(predict.ideology, probs = .05, na.rm = TRUE)
ci.upper.90 <- quantile(predict.ideology, probs = .95, na.rm = TRUE)

ideology.90 <- cbind(ci.low.90, mean, ci.upper.90)

plot.who<-data.frame(rbind(mid.95,up.95,col.95,ba.95, cambiemos.95,pj.95,
                           age.95,female.95,ideology.95,
                           mid.90,up.90,col.90,ba.90, cambiemos.90,pj.90,
                           age.90,female.90,ideology.90))


plot.who$variable<-c('Middle class','Upper class','College degree','Buenos Aires','Cambiemos identifier',
                     'PJ identifier','Age','Female','Ideology',
                     'Middle class','Upper class','College degree','Buenos Aires','Cambiemos identifier',
                     'PJ identifier','Age','Female','Ideology')

plot.who$ci<-c(rep('95',9),rep('90',9))



plot.who$variable <- factor(plot.who$variable, 
                            levels = c('Ideology','PJ identifier','Cambiemos identifier',
                                       'Age','Female',
                                       'College degree','Buenos Aires',
                                       'Upper class','Middle class'))


ggplot(plot.who,aes(x=variable,y=mean))+
  geom_hline(aes(yintercept=0), colour="#999999", linetype="dashed") +
  geom_pointrange( data =subset(plot.who,ci=='90'), stat = "identity", size =.6,position=pd,
                   mapping = aes(ymin= ci.low.95,ymax=ci.upper.95))+
  geom_pointrange( data =subset(plot.who,ci=='95'), stat = "identity", size =.3,position=pd,
                   mapping = aes(ymin= ci.low.95,ymax=ci.upper.95))+
  ylab("Change in probability of buying black market dollars") +xlab("Predictor") +
  theme_bw() +coord_flip()+
  #theme(aspect.ratio=.8)+
  theme(legend.position=c(.1,.9))+
  theme(axis.text.x  = element_text( size=10),axis.title.x = element_text(size=10)) +
  theme(axis.text.y  = element_text( size=10),axis.title.y = element_text(size=10))+
  theme(legend.background = element_rect( fill="white", size=.5, linetype="dotted"))+
  theme(legend.title = element_text(colour = 'black', size = 10))+
  theme(legend.text = element_text(colour = 'black', size = 8,vjust=0, hjust = 0))+
  theme(legend.key.width = unit(.5, "cm"))+
  theme(legend.key.size = unit(.5, "cm"))+
  theme(plot.margin = unit(c(.5, .5,.5, .5), "cm"))+
  
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
        panel.border = element_blank(),panel.background = element_blank())+
  theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())






# Appendix Table E3: Raw results from full MLE outcome model: Scioli vote

rows<-c('Intercept','Peronist Party ID','Cambiemos Party ID','Ideology','Middle class',
        'Upper class','Education','Female','Age','Amba')

tab.scioli<-data.frame(cbind(rows,scioli$par.treat,scioli$se.treat,scioli$par.control,scioli$se.control))
rows<-c(rows,'Sensitive item')
tab.scioli2<-data.frame(cbind(rows,scioli$par.outcome,scioli$se.outcome))
tab.scioli$id  <- 1:nrow(tab.scioli)
tab.scioli<-merge(tab.scioli,tab.scioli2,by.x='rows',by.y='rows',all.x=T,all.y=T)
tab.scioli<-subset(tab.scioli[order(tab.scioli$id), ],select=-(id))

colnames(tab.scioli)<-c('','Estimate','Std. error', 'Estimate','Std. error','Estimate','Std. error')
print(xtable(tab.scioli,digits=c(0,2,2,2,2,2,2,2)), include.rownames = F,sanitize.text.function=function(x){x})




# Table E4: Full MLE outcome model: Macri approval

macri<-  ictreg.joint(formula = list_outcome ~ pj+cambiemos+ideology+middleclass+upperclass+college+female+age+amba, data = df,
                      treat = "listgroup", outcome = "macri", J = 4, constrained = TRUE,
                      outcome.reg = "linear", maxIter = 1000)


rows<-c('Intercept','Peronist Party ID','Cambiemos Party ID','Ideology','Middle class',
        'Upper class','Education','Female','Age','Amba')

tab.macri<-data.frame(cbind(rows,macri$par.treat,macri$se.treat,macri$par.control,macri$se.control))
rows<-c(rows,'Sensitive item')
tab.macri2<-data.frame(cbind(rows,macri$par.outcome,macri$se.outcome))
tab.macri$id  <- 1:nrow(tab.macri)
tab.macri<-merge(tab.macri,tab.macri2,by.x='rows',by.y='rows',all.x=T,all.y=T)
tab.macri<-subset(tab.macri[order(tab.macri$id), ],select=-(id))

colnames(tab.macri)<-c('','Estimate','Std. error', 'Estimate','Std. error','Estimate','Std. error')
print(xtable(tab.macri,digits=c(0,2,2,2,2,2,2,2)), include.rownames = F,sanitize.text.function=function(x){x})





# Plots Figure 3. The effect of buying dollars on vote choice and presidential approval


# Predicted probabilities from outcome model: scioli vote
pred.scioli<-predict.ictreg.joint(scioli, se.fit = TRUE, interval = "confidence",
                                  level = 0.95, avg = TRUE,  constrained = TRUE,
                                  sensitive.value = "both",
                                  sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = TRUE)

x <- c(1,2,3,4)
avglist <- mean(c(pred.scioli$draws.mean[[1]], pred.scioli$draws.mean[[2]]))
avglist.lower <- quantile(c(pred.scioli$draws.mean[[1]], pred.scioli$draws.mean[[2]]), probs = .025)
avglist.upper <- quantile(c(pred.scioli$draws.mean[[1]], pred.scioli$draws.mean[[2]]), probs = .975)

#plot 
lwr <- c(avglist.lower, pred.scioli$fit[2, 2], pred.scioli$fit[1,2], pred.scioli$fit[3,2])
fit <- c(avglist, pred.scioli$fit[2, 1], pred.scioli$fit[1,1], pred.scioli$fit[3,1])
upr <- c(avglist.upper, pred.scioli$fit[2, 3], pred.scioli$fit[1,3], pred.scioli$fit[3,3])

data.scioli<-data.frame(cbind(fit,lwr,upr))
data.scioli$type<-factor(seq(1,4,1))

plot.scioli<-ggplot(subset(data.scioli,type!=4),aes(x=type,y=fit,colour=type,shape=type))+
  geom_point( data =subset(data.scioli,type!=4), stat = "identity", size =2,position=pd)+
  geom_errorbar( data =subset(data.scioli,type!=4), stat = "identity", size =.5,position=pd,width=0,
                 mapping = aes(ymin= lwr,ymax=upr))+
  
  ylab("Probability of voting for Scioli") +xlab("Participation in the black market") +
  theme_bw() +
  theme(legend.position='BLANK')+
  theme(axis.text.x  = element_text( size=10),axis.title.x = element_text(size=10)) +
  theme(axis.text.y  = element_text( size=10),axis.title.y = element_text(size=10))+
  theme(legend.background = element_rect( fill="white", size=.5, linetype="dotted"))+
  theme(legend.title = element_text(colour = 'black', size = 10))+
  theme(legend.text = element_text(colour = 'black', size = 8,vjust=0, hjust = 0))+
  theme(legend.key.width = unit(.5, "cm"))+
  theme(legend.key.size = unit(.5, "cm"))+
  theme(plot.margin = unit(c(.5, .5,.5, .5), "cm"))+
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank())+theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+
  scale_y_continuous(breaks=seq(0,.75,.25),labels=seq(0,.75,.25),lim=c(0,.75))+
  scale_x_discrete(breaks=seq(1,3,1),labels=c('Overall','Buyers','Non-buyers'))+
  annotate("text", x = 2, y = .7, label = paste('Difference Buyers and Non-buyers\n',
                                                round(data.scioli$fit[data.scioli$type==4],2),' [',
                                                round(data.scioli$lwr[data.scioli$type==4],2),',',
                                                round(data.scioli$upr[data.scioli$type==4],2),']'),size=3,colour='black')



# Predicted probabilities from outcome model: Macri approval

pred.macri<-predict.ictreg.joint(macri, se.fit = TRUE, interval = "confidence",
                                 level = 0.95, avg = TRUE,  constrained = TRUE,
                                 sensitive.value = "both",
                                 sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = TRUE)

x <- c(1,2,3,4)
avglist <- mean(c(pred.macri$draws.mean[[1]], pred.macri$draws.mean[[2]]))
avglist.lower <- quantile(c(pred.macri$draws.mean[[1]], pred.macri$draws.mean[[2]]), probs = .025)
avglist.upper <- quantile(c(pred.macri$draws.mean[[1]], pred.macri$draws.mean[[2]]), probs = .975)

#plot 
lwr <- c(avglist.lower, pred.macri$fit[2, 2], pred.macri$fit[1,2], pred.macri$fit[3,2])
fit <- c(avglist, pred.macri$fit[2, 1], pred.macri$fit[1,1], pred.macri$fit[3,1])
upr <- c(avglist.upper, pred.macri$fit[2, 3], pred.macri$fit[1,3], pred.macri$fit[3,3])

data.macri<-data.frame(cbind(fit,lwr,upr))

data.macri$type<-factor(seq(1,4,1))





plot.macri<-ggplot(subset(data.macri,type!=4),aes(x=type,y=fit,colour=type,shape=type))+
  geom_point( data =subset(data.macri,type!=4), stat = "identity", size =2,position=pd)+
  geom_errorbar( data =subset(data.macri,type!=4), stat = "identity", size =.5,position=pd,width=0,
                 mapping = aes(ymin= lwr,ymax=upr))+
  ylab("Presidential approval of Macri") +xlab("Participation in the black market") +
  theme_bw() +
  theme(legend.position='BLANK')+
  theme(axis.text.x  = element_text( size=10),axis.title.x = element_text(size=10)) +
  theme(axis.text.y  = element_text( size=10),axis.title.y = element_text(size=10))+
  theme(legend.background = element_rect( fill="white", size=.5, linetype="dotted"))+
  theme(legend.title = element_text(colour = 'black', size = 10))+
  theme(legend.text = element_text(colour = 'black', size = 8,vjust=0, hjust = 0))+
  theme(legend.key.width = unit(.5, "cm"))+
  theme(legend.key.size = unit(.5, "cm"))+
  theme(plot.margin = unit(c(.5, .5,.5, .5), "cm"))+
  
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank())+theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+
  scale_y_continuous(breaks=seq(1,4,1),labels=seq(1,4,1),lim=c(1,4))+
  scale_x_discrete(breaks=seq(1,3,1),labels=c('Overall','Buyers','Non-buyers'))+
  annotate("text", x = 2, y = 3.8, label = paste('Difference Buyers and Non-buyers\n',
                                                 round(data.macri$fit[data.macri$type==4],2),' [',
                                                 round(data.macri$lwr[data.macri$type==4],2),',',
                                                 round(data.macri$upr[data.macri$type==4],2),']'),size=3,colour='black')





grid.arrange(plot.scioli,plot.macri,nrow=1)

#### Figure 4: Effects of black market participation on knowledge about currency controls

# Support for cepo removal


cepo<-  ictreg.joint(formula = list_outcome ~ pj+cambiemos+ideology+middleclass+upperclass+college+female+age+amba, data = df,
                     treat = "listgroup", outcome = "cepo", J = 4, constrained = TRUE,
                     outcome.reg = "linear", maxIter = 1000)


# Knowledge: who adopted cepo


know_cepo<-  ictreg.joint(formula = list_outcome~ pj +cambiemos+ideology+middleclass+upperclass+college+female+age+amba, data = df,
                          treat = "listgroup", outcome = "know_cepo", J = 4, constrained = TRUE,
                          outcome.reg = "logistic", maxIter = 1000)

## Knowledge: who removed cepo



know_cepo2<-  ictreg.joint(formula = list_outcome~pj +cambiemos+ideology+middleclass+upperclass+college+female+age+amba, data = df,
                           treat = "listgroup", outcome = "know_cepo2", J = 4, constrained = TRUE,
                           outcome.reg = "logistic", maxIter = 1000)


################ Figure 4: The effect of buying dollars on policy preferences and knowledge################

# Predictred probabilities from outcome models
pred.cepo<-predict.ictreg.joint(cepo, se.fit = TRUE, interval = "confidence",
                                level = 0.95, avg = TRUE,  constrained = TRUE,
                                sensitive.value = "both",
                                sensitive.diff = TRUE, return.draws = TRUE,
                                predict.sensitive = TRUE)


pred.know.cepo<-predict.ictreg.joint(know_cepo, se.fit = TRUE, interval = "confidence",
                                     level = 0.95, avg = TRUE,  constrained = TRUE,
                                     sensitive.value = "both",
                                     sensitive.diff = TRUE, return.draws = TRUE,
                                     predict.sensitive = TRUE)


pred.know.cepo2<-predict.ictreg.joint(know_cepo2, se.fit = TRUE, interval = "confidence",
                                      level = 0.95, avg = TRUE,  constrained = TRUE,
                                      sensitive.value = "both",
                                      sensitive.diff = TRUE, return.draws = TRUE,
                                      predict.sensitive = TRUE)




psize<-2.5
textx<-10
texty<-10
titlex<-10
titley<-10
stext<-2.5
posy<-.99

x <- c(1,2,3,4)
avglist <- mean(c(pred.cepo$draws.mean[[1]], pred.cepo$draws.mean[[2]]))
avglist.lower <- quantile(c(pred.cepo$draws.mean[[1]], pred.cepo$draws.mean[[2]]), probs = .025)
avglist.upper <- quantile(c(pred.cepo$draws.mean[[1]], pred.cepo$draws.mean[[2]]), probs = .975)

#plot 
lwr <- c(avglist.lower, pred.cepo$fit[2, 2], pred.cepo$fit[1,2], pred.cepo$fit[3,2])
fit <- c(avglist, pred.cepo$fit[2, 1], pred.cepo$fit[1,1], pred.cepo$fit[3,1])
upr <- c(avglist.upper, pred.cepo$fit[2, 3], pred.cepo$fit[1,3], pred.cepo$fit[3,3])

data.cepo<-data.frame(cbind(fit,lwr,upr))

data.cepo$type<-factor(seq(1,4,1))

plot.cepo<-ggplot(subset(data.cepo,type!=4),aes(x=type,y=fit,colour=type,shape=type))+
  
  geom_point( data =subset(data.cepo,type!=4), stat = "identity", size=psize,position=pd)+
  geom_errorbar( data =subset(data.cepo,type!=4), stat = "identity", size=.5,position=pd,width=0,
                 mapping = aes(ymin= lwr,ymax=upr))+
  ylab("") +xlab("Participation in the black market") +
  theme_bw() + ggtitle("Support for cepo removal")+ 
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  theme(legend.position='BLANK')+
  theme(axis.text.x  = element_text( size=textx),axis.title.x = element_text(size=titlex)) +
  theme(axis.text.y  = element_text( size=texty),axis.title.y = element_text(size=titley))+
  theme(legend.background = element_rect( fill="white", size=.5, linetype="dotted"))+
  theme(legend.title = element_text(colour = 'black', size = 10))+
  theme(legend.text = element_text(colour = 'black', size = 8,vjust=0, hjust = 0))+
  theme(legend.key.width = unit(.5, "cm"))+
  theme(legend.key.size = unit(.5, "cm"))+
  theme(plot.margin = unit(c(.5, .5,.5, .5), "cm"))+
  
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank())+theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+
  scale_y_continuous(breaks=seq(1,5,1),labels=seq(1,5,1),lim=c(1,5))+
  scale_x_discrete(breaks=seq(1,3,1),labels=c('Overall','Buyers','Non-buyers'))+
  annotate("text", x = 2, y = 4.9, label = paste('Difference Buyers and Non-buyers\n',
                                                 round(data.cepo$fit[data.cepo$type==4],2),' [',
                                                 round(data.cepo$lwr[data.cepo$type==4],2),',',
                                                 round(data.cepo$upr[data.cepo$type==4],2),']')
           ,size=stext,colour='black')




avglist <- mean(c(pred.know.cepo$draws.mean[[1]], pred.know.cepo$draws.mean[[2]]))
avglist.lower <- quantile(c(pred.know.cepo$draws.mean[[1]], pred.know.cepo$draws.mean[[2]]), probs = .025)
avglist.upper <- quantile(c(pred.know.cepo$draws.mean[[1]], pred.know.cepo$draws.mean[[2]]), probs = .975)

#plot 
lwr <- c(avglist.lower, pred.know.cepo$fit[2, 2], pred.know.cepo$fit[1,2], pred.know.cepo$fit[3,2])
fit <- c(avglist, pred.know.cepo$fit[2, 1], pred.know.cepo$fit[1,1], pred.know.cepo$fit[3,1])
upr <- c(avglist.upper, pred.know.cepo$fit[2, 3], pred.know.cepo$fit[1,3], pred.know.cepo$fit[3,3])

data.know.cepo<-data.frame(cbind(fit,lwr,upr))

data.know.cepo$type<-factor(seq(1,4,1))


plot.know.cepo<-ggplot(subset(data.know.cepo,type!=4),aes(x=type,y=fit,colour=type,shape=type))+
  geom_point( data =subset(data.know.cepo,type!=4), stat = "identity", size=psize,position=pd)+
  geom_errorbar( data =subset(data.know.cepo,type!=4), stat = "identity", size=.5,position=pd,width=0,
                 mapping = aes(ymin= lwr,ymax=upr))+
  ylab("")+xlab("Participation in the black market") +
  theme_bw() + ggtitle("Knowledge: Who adopted cepo")+ 
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  theme(legend.position='BLANK')+
  theme(axis.text.x  = element_text( size=textx),axis.title.x = element_text(size=titlex)) +
  theme(axis.text.y  = element_text( size=texty),axis.title.y = element_text(size=titley))+
  theme(legend.background = element_rect( fill="white", size=.5, linetype="dotted"))+
  theme(legend.title = element_text(colour = 'black', size = 10))+
  theme(legend.text = element_text(colour = 'black', size = 8,vjust=0, hjust = 0))+
  theme(legend.key.width = unit(.5, "cm"))+
  theme(legend.key.size = unit(.5, "cm"))+
  theme(plot.margin = unit(c(.5, .5,.5, .5), "cm"))+
  
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank())+theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+
  scale_y_continuous(breaks=seq(0,1,.25),labels=seq(0,1,.25),lim=c(0,1))+
  scale_x_discrete(breaks=seq(1,3,1),labels=c('Overall','Buyers','Non-buyers'))+
  annotate("text", x = 2, y = posy, label = paste('Difference Buyers and Non-buyers\n',
                                                  round(data.know.cepo$fit[data.know.cepo$type==4],2),' [',
                                                  round(data.know.cepo$lwr[data.know.cepo$type==4],2),',',
                                                  round(data.know.cepo$upr[data.know.cepo$type==4],2),']')
           ,size=stext,colour='black')





avglist <- mean(c(pred.know.cepo2$draws.mean[[1]], pred.know.cepo2$draws.mean[[2]]))
avglist.lower <- quantile(c(pred.know.cepo2$draws.mean[[1]], pred.know.cepo2$draws.mean[[2]]), probs = .025)
avglist.upper <- quantile(c(pred.know.cepo2$draws.mean[[1]], pred.know.cepo2$draws.mean[[2]]), probs = .975)

#plot 
lwr <- c(avglist.lower, pred.know.cepo2$fit[2, 2], pred.know.cepo2$fit[1,2], pred.know.cepo2$fit[3,2])
fit <- c(avglist, pred.know.cepo2$fit[2, 1], pred.know.cepo2$fit[1,1], pred.know.cepo2$fit[3,1])
upr <- c(avglist.upper, pred.know.cepo2$fit[2, 3], pred.know.cepo2$fit[1,3], pred.know.cepo2$fit[3,3])

data.know.cepo2<-data.frame(cbind(fit,lwr,upr))

data.know.cepo2$type<-factor(seq(1,4,1))


plot.know.cepo2<-ggplot(subset(data.know.cepo2,type!=4),aes(x=type,y=fit,colour=type,shape=type))+
  geom_point( data =subset(data.know.cepo2,type!=4), stat = "identity", size=psize,position=pd)+
  geom_errorbar( data =subset(data.know.cepo2,type!=4), stat = "identity", size=.5,position=pd,width=0,
                 mapping = aes(ymin= lwr,ymax=upr))+
  theme_bw() + ggtitle("Knowledge: Who removed cepo")+ 
  ylab("") +xlab("Participation in the black market") +
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  theme(legend.position='BLANK')+
  theme(axis.text.x  = element_text( size=textx),axis.title.x = element_text(size=titlex)) +
  theme(axis.text.y  = element_text( size=texty),axis.title.y = element_text(size=titley))+
  theme(legend.background = element_rect( fill="white", size=.5, linetype="dotted"))+
  theme(legend.title = element_text(colour = 'black', size = 10))+
  theme(legend.text = element_text(colour = 'black', size = 8,vjust=0, hjust = 0))+
  theme(legend.key.width = unit(.5, "cm"))+
  theme(legend.key.size = unit(.5, "cm"))+
  theme(plot.margin = unit(c(.5, .5,.5, .5), "cm"))+
  
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank())+theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+
  scale_y_continuous(breaks=seq(0,1,.25),labels=seq(0,1,.25),lim=c(0,1))+
  scale_x_discrete(breaks=seq(1,3,1),labels=c('Overall','Buyers','Non-buyers'))+
  annotate("text", x = 2, y = posy, label = paste('Difference Buyers and Non-buyers\n',
                                                  round(data.know.cepo2$fit[data.know.cepo2$type==4],2),' [',
                                                  round(data.know.cepo2$lwr[data.know.cepo2$type==4],2),',',
                                                  round(data.know.cepo2$upr[data.know.cepo2$type==4],2),']'),size=stext,colour='black')


grid.arrange(plot.cepo,plot.know.cepo,plot.know.cepo2,nrow=1)


#######  Raw results from models with knowledge and policy preference outcomes ####### 

# Appendix Table E5
rows<-c('Intercept','Peronist Party ID','Cambiemos Party ID','Ideology','Middle class',
        'Upper class','Education','Female','Age','Amba')

tab.cepo<-data.frame(cbind(rows,cepo$par.treat,cepo$se.treat,cepo$par.control,cepo$se.control))
rows<-c(rows,'Sensitive item')
tab.cepo2<-data.frame(cbind(rows,cepo$par.outcome,cepo$se.outcome))
tab.cepo$id  <- 1:nrow(tab.cepo)
tab.cepo<-merge(tab.cepo,tab.cepo2,by.x='rows',by.y='rows',all.x=T,all.y=T)
tab.cepo<-subset(tab.cepo[order(tab.cepo$id), ],select=-(id))

colnames(tab.cepo)<-c('','Estimate','Std. error', 'Estimate','Std. error','Estimate','Std. error')
print(xtable(tab.cepo,digits=c(0,2,2,2,2,2,2,2)), include.rownames = F,sanitize.text.function=function(x){x})


# Appendix Table E6
rows<-c('Intercept','Peronist Party ID','Cambiemos Party ID','Ideology','Middle class',
        'Upper class','Education','Female','Age','Amba')

tab.know_cepo<-data.frame(cbind(rows,know_cepo$par.treat,know_cepo$se.treat,know_cepo$par.control,know_cepo$se.control))
rows<-c(rows,'Sensitive item')
tab.know_cepo2<-data.frame(cbind(rows,know_cepo$par.outcome,know_cepo$se.outcome))
tab.know_cepo$id  <- 1:nrow(tab.know_cepo)
tab.know_cepo<-merge(tab.know_cepo,tab.know_cepo2,by.x='rows',by.y='rows',all.x=T,all.y=T)
tab.know_cepo<-subset(tab.know_cepo[order(tab.know_cepo$id), ],select=-(id))

colnames(tab.know_cepo)<-c('','Estimate','Std. error', 'Estimate','Std. error','Estimate','Std. error')
print(xtable(tab.know_cepo,digits=c(0,2,2,2,2,2,2,2)), include.rownames = F,sanitize.text.function=function(x){x})


# Appendix Table E7

rows<-c('Intercept','Peronist Party ID','Cambiemos Party ID','Ideology','Middle class',
        'Upper class','Education','Female','Age','Amba')

tab.know_cepo2<-data.frame(cbind(rows,know_cepo2$par.treat,know_cepo2$se.treat,know_cepo2$par.control,know_cepo2$se.control))
rows<-c(rows,'Sensitive item')
tab.know_cepo22<-data.frame(cbind(rows,know_cepo2$par.outcome,know_cepo2$se.outcome))
tab.know_cepo2$id  <- 1:nrow(tab.know_cepo2)
tab.know_cepo2<-merge(tab.know_cepo2,tab.know_cepo22,by.x='rows',by.y='rows',all.x=T,all.y=T)
tab.know_cepo2<-subset(tab.know_cepo2[order(tab.know_cepo2$id), ],select=-(id))

colnames(tab.know_cepo2)<-c('','Estimate','Std. error', 'Estimate','Std. error','Estimate','Std. error')
print(xtable(tab.know_cepo2,digits=c(0,2,2,2,2,2,2,2)), include.rownames = F,sanitize.text.function=function(x){x})



# Placebo outcomes

# Knowledge about CCTs

know_auh<-  ictreg.joint(formula = list_outcome~pj +cambiemos+ideology+middleclass+upperclass+college+female+age+amba, data = df,
                         treat = "listgroup", outcome = "know_auh", J = 4, constrained = TRUE,
                         outcome.reg = "logistic", maxIter = 1000)


# Appendix Table E8

rows<-c('Intercept','Peronist Party ID','Cambiemos Party ID','Ideology','Middle class',
        'Upper class','Education','Female','Age','Amba')

tab.know_auh<-data.frame(cbind(rows,know_auh$par.treat,know_auh$se.treat,know_auh$par.control,know_auh$se.control))
rows<-c(rows,'Sensitive item')
tab.know_auh2<-data.frame(cbind(rows,know_auh$par.outcome,know_auh$se.outcome))
tab.know_auh$id  <- 1:nrow(tab.know_auh)
tab.know_auh<-merge(tab.know_auh,tab.know_auh2,by.x='rows',by.y='rows',all.x=T,all.y=T)
tab.know_auh<-subset(tab.know_auh[order(tab.know_auh$id), ],select=-(id))

colnames(tab.know_auh)<-c('','Estimate','Std. error', 'Estimate','Std. error','Estimate','Std. error')
print(xtable(tab.know_auh,digits=c(0,2,2,2,2,2,2,2)), include.rownames = F,sanitize.text.function=function(x){x})


# Knowledge about the nationalization of YPF
know_ypf<-  ictreg.joint(formula = list_outcome~ pj+cambiemos+ideology+middleclass+upperclass+college+female+age+amba, data = df,
                         treat = "listgroup", outcome = "know_ypf", J = 4, constrained = TRUE,
                         outcome.reg = "logistic", maxIter = 1000)

rows<-c('Intercept','Peronist Party ID','Cambiemos Party ID','Ideology','Middle class',
        'Upper class','Education','Female','Age','Amba')

# Appendix Table E9
tab.know_ypf<-data.frame(cbind(rows,know_ypf$par.treat,know_ypf$se.treat,know_ypf$par.control,know_ypf$se.control))
rows<-c(rows,'Sensitive item')
tab.know_ypf2<-data.frame(cbind(rows,know_ypf$par.outcome,know_ypf$se.outcome))
tab.know_ypf$id  <- 1:nrow(tab.know_ypf)
tab.know_ypf<-merge(tab.know_ypf,tab.know_ypf2,by.x='rows',by.y='rows',all.x=T,all.y=T)
tab.know_ypf<-subset(tab.know_ypf[order(tab.know_ypf$id), ],select=-(id))

colnames(tab.know_ypf)<-c('','Estimate','Std. error', 'Estimate','Std. error','Estimate','Std. error')
print(xtable(tab.know_ypf,digits=c(0,2,2,2,2,2,2,2)), include.rownames = F,sanitize.text.function=function(x){x})


# Knowledge about the adoption of gay marriage law
know_gay<-  ictreg.joint(formula = list_outcome~pj +cambiemos+ideology+middleclass+upperclass+college+female+age+amba, data = df,
                         treat = "listgroup", outcome = "know_gay", J = 4, constrained = TRUE,
                         outcome.reg = "logistic", maxIter = 1000)


# Appendix Table E10
rows<-c('Intercept','Peronist Party ID','Cambiemos Party ID','Ideology','Middle class',
        'Upper class','Education','Female','Age','Amba')

tab.know_gay<-data.frame(cbind(rows,know_gay$par.treat,know_gay$se.treat,know_gay$par.control,know_gay$se.control))
rows<-c(rows,'Sensitive item')
tab.know_gay2<-data.frame(cbind(rows,know_gay$par.outcome,know_gay$se.outcome))
tab.know_gay$id  <- 1:nrow(tab.know_gay)
tab.know_gay<-merge(tab.know_gay,tab.know_gay2,by.x='rows',by.y='rows',all.x=T,all.y=T)
tab.know_gay<-subset(tab.know_gay[order(tab.know_gay$id), ],select=-(id))

colnames(tab.know_gay)<-c('','Estimate','Std. error', 'Estimate','Std. error','Estimate','Std. error')
print(xtable(tab.know_gay,digits=c(0,2,2,2,2,2,2,2)), include.rownames = F,sanitize.text.function=function(x){x})






############ Figure E4. Placebo tests #################

# Predicted probabilities from models with placebo outcomes

pred.know.auh<-predict.ictreg.joint(know_auh, se.fit = TRUE, interval = "confidence",
                                    level = 0.95, avg = TRUE,  constrained = TRUE,
                                    sensitive.value = "both",
                                    sensitive.diff = TRUE, return.draws = TRUE,
                                    predict.sensitive = TRUE)


pred.know.ypf<-predict.ictreg.joint(know_ypf, se.fit = TRUE, interval = "confidence",
                                    level = 0.95, avg = TRUE,  constrained = TRUE,
                                    sensitive.value = "both",
                                    sensitive.diff = TRUE, return.draws = TRUE,
                                    predict.sensitive = TRUE)



pred.know.gay<-predict.ictreg.joint(know_gay, se.fit = TRUE, interval = "confidence",
                                    level = 0.95, avg = TRUE,  constrained = TRUE,
                                    sensitive.value = "both",
                                    sensitive.diff = TRUE, return.draws = TRUE,
                                    predict.sensitive = TRUE)



avglist <- mean(c(pred.know.ypf$draws.mean[[1]], pred.know.ypf$draws.mean[[2]]))
avglist.lower <- quantile(c(pred.know.ypf$draws.mean[[1]], pred.know.ypf$draws.mean[[2]]), probs = .025)
avglist.upper <- quantile(c(pred.know.ypf$draws.mean[[1]], pred.know.ypf$draws.mean[[2]]), probs = .975)

#plot 
lwr <- c(avglist.lower, pred.know.ypf$fit[2, 2], pred.know.ypf$fit[1,2], pred.know.ypf$fit[3,2])
fit <- c(avglist, pred.know.ypf$fit[2, 1], pred.know.ypf$fit[1,1], pred.know.ypf$fit[3,1])
upr <- c(avglist.upper, pred.know.ypf$fit[2, 3], pred.know.ypf$fit[1,3], pred.know.ypf$fit[3,3])

data.know.ypf<-data.frame(cbind(fit,lwr,upr))

data.know.ypf$type<-factor(seq(1,4,1))


plot.know.ypf<-ggplot(subset(data.know.ypf,type!=4),aes(x=type,y=fit,colour=type,shape=type))+
  
  geom_pointrange( data =subset(data.know.ypf,type!=4), stat = "identity", size=.8,position=pd,
                   mapping = aes(ymin= lwr,ymax=upr))+
  theme_bw() + ggtitle("Which president nationalized oil company?")+ 
  ylab("") +xlab("Participation in the black market") +
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  theme(legend.position='BLANK')+
  theme(axis.text.x  = element_text( size=12),axis.title.x = element_text(size=12)) +
  theme(axis.text.y  = element_text( size=8),axis.title.y = element_text(size=12))+
  theme(legend.background = element_rect( fill="white", size=.5, linetype="dotted"))+
  theme(legend.title = element_text(colour = 'black', size = 10))+
  theme(legend.text = element_text(colour = 'black', size = 8,vjust=0, hjust = 0))+
  theme(legend.key.width = unit(.5, "cm"))+
  theme(legend.key.size = unit(.5, "cm"))+
  theme(plot.margin = unit(c(.5, .5,.5, .5), "cm"))+
  
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank())+theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+
  scale_y_continuous(breaks=seq(0,1,.25),labels=seq(0,1,.25),lim=c(0,1))+
  scale_x_discrete(breaks=seq(1,3,1),labels=c('Overall','Buyers','Non-buyers'))+
  annotate("text", x = 2, y = .98, label = paste('Diff Buyers Non-buyers\n',
                                                 round(data.know.ypf$fit[data.know.ypf$type==4],2),' [',
                                                 round(data.know.ypf$lwr[data.know.ypf$type==4],2),',',
                                                 round(data.know.ypf$upr[data.know.ypf$type==4],2),']')
           ,size=3,colour='black')

avglist <- mean(c(pred.know.auh$draws.mean[[1]], pred.know.auh$draws.mean[[2]]))
avglist.lower <- quantile(c(pred.know.auh$draws.mean[[1]], pred.know.auh$draws.mean[[2]]), probs = .025)
avglist.upper <- quantile(c(pred.know.auh$draws.mean[[1]], pred.know.auh$draws.mean[[2]]), probs = .975)

#plot 
lwr <- c(avglist.lower, pred.know.auh$fit[2, 2], pred.know.auh$fit[1,2], pred.know.auh$fit[3,2])
fit <- c(avglist, pred.know.auh$fit[2, 1], pred.know.auh$fit[1,1], pred.know.auh$fit[3,1])
upr <- c(avglist.upper, pred.know.auh$fit[2, 3], pred.know.auh$fit[1,3], pred.know.auh$fit[3,3])

data.know.auh<-data.frame(cbind(fit,lwr,upr))

data.know.auh$type<-factor(seq(1,4,1))


plot.know.auh<-ggplot(subset(data.know.auh,type!=4),aes(x=type,y=fit,colour=type,shape=type))+
  
  geom_pointrange( data =subset(data.know.auh,type!=4), stat = "identity", size=.8,position=pd,
                   mapping = aes(ymin= lwr,ymax=upr))+
  theme_bw() + ggtitle("Which president adopted CCT program?")+ 
  ylab("") +xlab("Participation in the black market") +
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  theme(legend.position='BLANK')+
  theme(axis.text.x  = element_text( size=12),axis.title.x = element_text(size=12)) +
  theme(axis.text.y  = element_text( size=8),axis.title.y = element_text(size=12))+
  theme(legend.background = element_rect( fill="white", size=.5, linetype="dotted"))+
  theme(legend.title = element_text(colour = 'black', size = 10))+
  theme(legend.text = element_text(colour = 'black', size = 8,vjust=0, hjust = 0))+
  theme(legend.key.width = unit(.5, "cm"))+
  theme(legend.key.size = unit(.5, "cm"))+
  theme(plot.margin = unit(c(.5, .5,.5, .5), "cm"))+
  
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank())+theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+
  scale_y_continuous(breaks=seq(0,1,.25),labels=seq(0,1,.25),lim=c(0,1))+
  scale_x_discrete(breaks=seq(1,3,1),labels=c('Overall','Buyers','Non-buyers'))+
  annotate("text", x = 2, y = .98, label = paste('Diff Buyers Non-buyers\n',
                                                 round(data.know.auh$fit[data.know.auh$type==4],2),' [',
                                                 round(data.know.auh$lwr[data.know.auh$type==4],2),',',
                                                 round(data.know.auh$upr[data.know.auh$type==4],2),']')
           ,size=3,colour='black')


avglist <- mean(c(pred.know.merc$draws.mean[[1]], pred.know.merc$draws.mean[[2]]))
avglist.lower <- quantile(c(pred.know.merc$draws.mean[[1]], pred.know.merc$draws.mean[[2]]), probs = .025)
avglist.upper <- quantile(c(pred.know.merc$draws.mean[[1]], pred.know.merc$draws.mean[[2]]), probs = .975)

#plot 
lwr <- c(avglist.lower, pred.know.merc$fit[2, 2], pred.know.merc$fit[1,2], pred.know.merc$fit[3,2])
fit <- c(avglist, pred.know.merc$fit[2, 1], pred.know.merc$fit[1,1], pred.know.merc$fit[3,1])
upr <- c(avglist.upper, pred.know.merc$fit[2, 3], pred.know.merc$fit[1,3], pred.know.merc$fit[3,3])


avglist <- mean(c(pred.know.gay$draws.mean[[1]], pred.know.gay$draws.mean[[2]]))
avglist.lower <- quantile(c(pred.know.gay$draws.mean[[1]], pred.know.gay$draws.mean[[2]]), probs = .025)
avglist.upper <- quantile(c(pred.know.gay$draws.mean[[1]], pred.know.gay$draws.mean[[2]]), probs = .975)

#plot 
lwr <- c(avglist.lower, pred.know.gay$fit[2, 2], pred.know.gay$fit[1,2], pred.know.gay$fit[3,2])
fit <- c(avglist, pred.know.gay$fit[2, 1], pred.know.gay$fit[1,1], pred.know.gay$fit[3,1])
upr <- c(avglist.upper, pred.know.gay$fit[2, 3], pred.know.gay$fit[1,3], pred.know.gay$fit[3,3])

data.know.gay<-data.frame(cbind(fit,lwr,upr))

data.know.gay$type<-factor(seq(1,4,1))


plot.know.gay<-ggplot(subset(data.know.gay,type!=4),aes(x=type,y=fit,colour=type,shape=type))+
  
  geom_pointrange( data =subset(data.know.gay,type!=4), stat = "identity", size=.8,position=pd,
                   mapping = aes(ymin= lwr,ymax=upr))+
  ylab("") +xlab("") +
  theme_bw() + ggtitle("Which president adopted gay marriage law?")+ 
  ylab("") +xlab("Participation in the black market") +
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  #theme(aspect.ratio=.8)+
  theme(legend.position='BLANK')+
  theme(axis.text.x  = element_text( size=12),axis.title.x = element_text(size=12)) +
  theme(axis.text.y  = element_text( size=8),axis.title.y = element_text(size=12))+
  theme(legend.background = element_rect( fill="white", size=.5, linetype="dotted"))+
  theme(legend.title = element_text(colour = 'black', size = 10))+
  theme(legend.text = element_text(colour = 'black', size = 8,vjust=0, hjust = 0))+
  theme(legend.key.width = unit(.5, "cm"))+
  theme(legend.key.size = unit(.5, "cm"))+
  theme(plot.margin = unit(c(.5, .5,.5, .5), "cm"))+
  
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank())+theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+
  scale_y_continuous(breaks=seq(0,1,.25),labels=seq(0,1,.25),lim=c(0,1))+
  scale_x_discrete(breaks=seq(1,3,1),labels=c('Overall','Buyers','Non-buyers'))+
  annotate("text", x = 2, y = .98, label = paste('Diff Buyers Non-buyers\n',
                                                 round(data.know.gay$fit[data.know.gay$type==4],2),' [',
                                                 round(data.know.gay$lwr[data.know.gay$type==4],2),',',
                                                 round(data.know.gay$upr[data.know.gay$type==4],2),']'),size=3,colour='black')



grid.arrange(plot.know.ypf,plot.know.auh,plot.know.gay,nrow=1)


###### Figure E5. Comparing Effects of Sensitive item to other covariates on vote for Scioli


vscioli.sens <- pred.scioli$fit[3, ]


#### Now we estimate difference in predicted probability of voting for Scioli for
### different values of covariates included in the outcome model

# Education

vscioli.edu<-predict.ictreg.joint(scioli,newdata=subset(df,college==0),newdata.diff=subset(df,college==1),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value =0, sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = F)


pred.vscioli.edu <- vscioli.edu$draws.mean[[2]] -  vscioli.edu$draws.mean[[1]]
pred.vscioli.edu <- c(mean(pred.vscioli.edu),quantile(pred.vscioli.edu, probs = .025, na.rm = TRUE),
                      quantile(pred.vscioli.edu, probs = .975, na.rm = TRUE))

# Middle class

vscioli.mid<-predict.ictreg.joint(scioli,newdata=subset(df,middleclass==0),newdata.diff=subset(df,middleclass==1),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value =0, sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = F)


pred.vscioli.mid <- vscioli.mid$draws.mean[[2]] -  vscioli.mid$draws.mean[[1]]
pred.vscioli.mid <- c(mean(pred.vscioli.mid),quantile(pred.vscioli.mid, probs = .025, na.rm = TRUE),
                      quantile(pred.vscioli.mid, probs = .975, na.rm = TRUE))


# Upper class
vscioli.upp<-predict.ictreg.joint(scioli,newdata=subset(df,upperclass==0),newdata.diff=subset(df,upperclass==1),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value =0, sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = F)


pred.vscioli.upp <- vscioli.upp$draws.mean[[2]] -  vscioli.upp$draws.mean[[1]]
pred.vscioli.upp <- c(mean(pred.vscioli.upp),quantile(pred.vscioli.upp, probs = .025, na.rm = TRUE),
                      quantile(pred.vscioli.upp, probs = .975, na.rm = TRUE))

# Cambemos party ID

vscioli.camb<-predict.ictreg.joint(scioli,newdata=subset(df,cambiemos==0),newdata.diff=subset(df,cambiemos==1),
                                   se.fit = TRUE, interval = "confidence", level = 0.95,
                                   avg = TRUE, sensitive.value =0, sensitive.diff = TRUE, return.draws = TRUE,
                                   predict.sensitive = F)


pred.vscioli.camb <- vscioli.camb$draws.mean[[2]] -  vscioli.camb$draws.mean[[1]]
pred.vscioli.camb <- c(mean(pred.vscioli.camb),quantile(pred.vscioli.camb, probs = .025, na.rm = TRUE),
                       quantile(pred.vscioli.camb, probs = .975, na.rm = TRUE))


# PJ party ID
vscioli.pj<-predict.ictreg.joint(scioli,newdata=subset(df,pj==0),newdata.diff=subset(df,pj==1),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value =0, sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = F)


pred.vscioli.pj <- vscioli.pj$draws.mean[[2]] -  vscioli.pj$draws.mean[[1]]
pred.vscioli.pj <- c(mean(pred.vscioli.pj),quantile(pred.vscioli.pj, probs = .025, na.rm = TRUE),
                     quantile(pred.vscioli.pj, probs = .975, na.rm = TRUE))



# Ideology
vscioli.id<-predict.ictreg.joint(scioli,newdata=subset(df,ideology=5),newdata.diff=subset(df,ideology==7),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value =0, sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = F)


pred.vscioli.id <- vscioli.id$draws.mean[[2]] -  vscioli.id$draws.mean[[1]]
pred.vscioli.id <- c(mean(pred.vscioli.id),quantile(pred.vscioli.id, probs = .025, na.rm = TRUE),
                     quantile(pred.vscioli.id, probs = .975, na.rm = TRUE))


# Female
vscioli.fem<-predict.ictreg.joint(scioli,newdata=subset(df,female==0),newdata.diff=subset(df,female==1),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value =0, sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = F)


pred.vscioli.fem <- vscioli.fem$draws.mean[[2]] -  vscioli.fem$draws.mean[[1]]
pred.vscioli.fem <- c(mean(pred.vscioli.fem),quantile(pred.vscioli.fem, probs = .025, na.rm = TRUE),
                      quantile(pred.vscioli.fem, probs = .975, na.rm = TRUE))


# Age
vscioli.age<-predict.ictreg.joint(scioli,newdata=subset(df,age==1),newdata.diff=subset(df,age==3),
                                  se.fit = TRUE, interval = "confidence", level = 0.95,
                                  avg = TRUE, sensitive.value =0, sensitive.diff = TRUE, return.draws = TRUE,
                                  predict.sensitive = F)


pred.vscioli.age <- vscioli.age$draws.mean[[2]] -  vscioli.age$draws.mean[[1]]
pred.vscioli.age <- c(mean(pred.vscioli.age),quantile(pred.vscioli.age, probs = .025, na.rm = TRUE),
                      quantile(pred.vscioli.age, probs = .975, na.rm = TRUE))


# BA region
vscioli.ba<-predict.ictreg.joint(scioli,newdata=subset(df,amba==0),newdata.diff=subset(df,amba==1),
                                 se.fit = TRUE, interval = "confidence", level = 0.95,
                                 avg = TRUE, sensitive.value =0, sensitive.diff = TRUE, return.draws = TRUE,
                                 predict.sensitive = F)


pred.vscioli.ba <- vscioli.ba$draws.mean[[2]] -  vscioli.ba$draws.mean[[1]]
pred.vscioli.ba <- c(mean(pred.vscioli.ba),quantile(pred.vscioli.ba, probs = .025, na.rm = TRUE),
                     quantile(pred.vscioli.ba, probs = .975, na.rm = TRUE))








plot.size<-data.frame(rbind(vscioli.sens,pred.vscioli.edu,pred.vscioli.mid,pred.vscioli.upp,
                            pred.vscioli.camb,pred.vscioli.pj,pred.vscioli.id,
                            pred.vscioli.ba,pred.vscioli.age,pred.vscioli.fem))


plot.size$variable<-c('Sensitive item','College degree','Middle class','Upper class',
                      'Cambiemos identifier','PJ identifier','Ideology','Buenos Aires','Age','Female')


plot.size$variable <- factor(plot.size$variable, 
                             levels = c('Sensitive item','Ideology','PJ identifier','Cambiemos identifier',
                                        'Age','Female',
                                        'College degree','Buenos Aires',
                                        'Upper class','Middle class'))

plot.size
pd <- position_dodge(.5)

ggplot(plot.size,aes(x=variable,y=fit))+
  geom_hline(aes(yintercept=0), colour="#999999", linetype="dashed") +
  geom_pointrange( data =plot.size, stat = "identity", size =.6,position=pd,
                   mapping = aes(ymin= lwr,ymax=upr))+
  ylab("Change in probability of voting for Scioli") +xlab("Predictor") +
  theme_bw() +coord_flip()+
  #theme(aspect.ratio=.8)+
  theme(legend.position=c(.1,.9))+
  theme(axis.text.x  = element_text( size=10),axis.title.x = element_text(size=10)) +
  theme(axis.text.y  = element_text( size=10),axis.title.y = element_text(size=10))+
  theme(legend.background = element_rect( fill="white", size=.5, linetype="dotted"))+
  theme(legend.title = element_text(colour = 'black', size = 10))+
  theme(legend.text = element_text(colour = 'black', size = 8,vjust=0, hjust = 0))+
  theme(legend.key.width = unit(.5, "cm"))+
  theme(legend.key.size = unit(.5, "cm"))+
  theme(plot.margin = unit(c(.5, .5,.5, .5), "cm"))+
  
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
        panel.border = element_blank(),panel.background = element_blank())+
  theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())



